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SECTION  1 
INTRODUCTION 

The  prediction  of  blast  wave  environments,  whether  a  result  of  a  chemical  or  a  nuclear 
source,  is  an  extremely  challenging  technical  problem.  A  combination  of  analytical  reason¬ 
ing,  study  of  the  records  of  nonrepeatable  historical  events,  design  and  study  of  (perhaps) 
repeatable  large-scale  field  simulations,  laboratory  experiments  and  studies  involving  com¬ 
putational  fluid  dynamics  (CFD)  are  all  involved  in  a  highly  interactive  and  fast-changing 
dynamic  process.  Among  the  difficulties  encountered  are  (1)  the  problems  posed  are  often 
three-dimensional  whereas  simulations  are  more  easily  arranged  for  two-dimensional  anal¬ 
ysis,  (2)  many  unusual  physical  phenomena  must  be  dealt  with  including,  but  not  limited 
to,  high  temperature  real  gas  effects,  very  strong  incident  shock  waves,  rapid  fireball  rise 
leading  to  a  multiscale  flowfield,  complex  and  unusual  multiphase  boundary  layer  effects 
including  cratering  and  precursor/decursor  effects  in  certain  situations  and  (3)  the  neces¬ 
sity  of  multimaterial  analysis  in  certain  situations,  e.g.,  the  air/product  gas  interface  for  a 
chemical  explosive  or  problems  involving  water. 

This  report  is  meant  to  summarize  the  progress  made  during  the  period  of  (approx¬ 
imately)  1987-1990  in  developing  a  particular  CFD  technology  based  on  a  second  order 
adaptation  of  the  ideas  of  Godunov  [33],  [34];  also,  we  discuss  the  many  applications  that 
have  been  studied  using  these  new  methods  during  this  period.  The  current  state-of-the- 
art  is  assessed  as  well. 

In  order  to  set  the  stage  for  such  a  discussion,  we  first  summarize  the  state  of  affairs 
in  the  relevant  CFD  community  at  around  1978  when  I  and  others  involved  in  the  work 
here  began  to  look  at  a  certain  set  of  experimental  data  which  did  not  seem  amenable  to 
CFD  analysis.  At  that  time,  all  of  the  heavily  used  CFD  codes  were  based  on  the  ideas  of 
the  Lax-Wendroff  scheme.  The  aerodynamic  community  used  a  relatively  clean  version  - 
MacCormack’s  scheme  -  and  were  mainly  interested  in  developing,  e.g.,  implicit  versions  to 
quickly  attain  steady  state  flows.  The  Defense  Nuclear  Agency  (DNA)  community,  on  the 
other  hand,  used  codes  such  as  the  HULL  code  which  contained  a  substantially  enhanced 
set  of  capabilities  to  enable  it  to  handle  the  range  of  phenomonology  discussed  above; 
of  course,  steady  state  convergence  was  not  an  issue.  Additionally,  many  other  codes 
developed  at  the  weapons  labs  (Los  Alamos  National  Laboratory  (LANL)  and  Lawrence 
Livermore  National  Laboratory  (LLNL),  today)  were  also  used.  For  example,  nuclear 
weapon  design  obviously  requires  multimaterial  capability,  and  this  was  available  since  at 
least  the  early  1960’s.  During  this  period,  the  first  really  new  method  -  the  Flux  Corrected 
Transport  (FCT)  scheme  developed  at  Naval  Research  Laboratory  (NRL)  -  was  also  being 
used  for  unsteady  work.  A  very  important  observation  is  that  the  speed  and  memory  of 
the  pre-CRAY  machines  were  a  tiny  fraction  of  the  capabilities  that  we  take  for  granted 
today;  the  same  can  be  said  for  ease  of  use,  post-processing,  etc.  As  a  result,  it  was 
not  possible  to  make  truly  high  resolution  (i.e.,  in  some  sense,  converged)  calculations  for 
difficult  problems;  hence,  it  really  was  not  known  whether  or  not  the  information  obtained 


1 


from  them  was  reliable.  So,  just  as  the  aerodynamic  community  relied  mainly  on  wind 
tunnels,  the  DNA  community  used  field  tests. 

At  this  time,  a  new  problem  arose  in  the  field  of  high  overpressure  shock  and  b<^t 
waves.  To  be  precise,  experimentalists  in  both  the  laboratory  and  the  field  were  observing, 
for  certain  cases,  multiple  pressure  peaks  behind  the  incident  wave;  the  CFD  codes  were 
unable  to  predict  the  phenomenon  and  there  did  not  seem  to  be  a  viable  analytical  theory 
to  explain  the  extra  peaks.  Simultaneously,  three  events  took  place  which  were  to  resolve 
the  problem  about  three  years  later:  (1)  the  CRAY  1  computer  became  available,  (i)  Prof. 
I.  I.  Glass  at  University  of  Toronto  Institute  of  Aerospace  Studies  (UTIAS)  began  using 
infinite-fringe  interferometry  in  analyzing  self-similar  laboratory  scale  shock  wave  inter¬ 
actions  and  (3)  the  second  order  Godunov  scheme  was  invented  by  van  Leer  [64]  (and  the 
FCT  group  was  also  making  rapid  progress  with  related  ideas).  By  1982,  it  was  completely 
clear  that,  at  least  for  the  laboratory  scale  experiments,  the  extra  peak(s)  were  real  and 
were  due  to  the  complexities  of  double  Mach  stems.  This  was  independently  discovered  by 
the  UTIAS  interferograms,  the  FCT  calculations  and  the  second  order  Godunov  calcula¬ 
tions.  At  this  time,  the  UTIAS  group  began  a  lengthy  collaboration  with  the  Naval  Surface 
Warfare  Center  (NSWC)  /  Lawrence  Berkeley  Laboratory  (LBL)  group  which  culminated 
in  a  direct  comparison  study  (laboratory  experiments  and  calculations)  using  real  (this  is 
obvious  for  the  experiments  but  difficult  for  computations  because  nontrivial  equation-of- 
state  (EOS)  software  must  be  used)  air;  the  extensive  agreement  between  the  two  sets  of 
results  served  to  confirm  both  groups’  methods.  The  computations  were  performed  with 
a  variant  of  van  Leer’s  MUSCL  scheme  due  to  Colella  and  Woodward  [16],  [69]  as  further 
modified  by  Colella  and  Glaz  [13];  the  latter  code  was  the  first  with  general  (convex  only, 
however)  EOS  capability.  The  results  were  published  a  few  years  later,  see  [28],  [29].  At 
this  point,  the  FCT  (and  even  the  HULL)  results  were  also  confirming  the  basic  ideas.  A 
bit  later,  another  comparison  study  was  undertaken  using  SF6  as  the  test  gas,  see  [30]. 
Since  the  assumption  of  equilibrium  is  much  closer  to  being  correct  for  SFg  than  for  low 
density  air,  these  comparisons  were  even  better;  indeed,  they  essentially  match  exactly.  In 
Section  1  below,  some  further  details  and  new  work  are  briefly  discussed  for  the  self-similar 
case. 

Simultaneously  with  this  work  aimed  at  the  self-similar  case,  many  groups  were  now 
going  after  high  resolution  results  for  height-of-burst  studies  which  are  truly  unsteady 
and  considerably  more  computationally  expensive.  The  NSWC/LBL  group  chose  a  set 
of  experiments  which  were  carefully  analyzed  by  J.  Carpenter.  The  calculations,  which 
are  started  by  fully  burning  the  appropriate  high  explosive  in  an  axisymmetric  geometry, 
were  performed  and  analyzed  by  Colella,  et  al.  [12].  The  results  matched  extremely  well 
with  the  available  data  and,  more  importantly,  verified  the  multiple  peak  phenomenon;  in 
fact,  three  distinct  peaks  at  close-in  stations  were  observed  in  the  most  highly  resolved 
computation.  This,  and  later  CFD  work  by  various  groups  finally  settled  the  original 
problem. 
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At  this  point,  a  new  set  of  problems  became  of  interest.  An  intriguing  set  of  pho¬ 
tographs  from  early  nuclear  e”^'.  showed  a  clear  precursor  wave  moving  well  out  in  front 
of  the  incident  blast  wave.  Here,  an  analytic  hypothesis  could  and  was  formed  to  explain 
the  phenomenon:  radiation  heating  of  the  ground  and  nearby  boundary  layer  air  would 
increase  the  shock  speed  near  the  ground  but  not  too  far  above  it.  Also,  applications  were 
leading  the  CFD  groups  into  the  study  of  blast  wave  phenomena  with  nonideal  boundary 
conditions  and  various  modeling  hypotheses  were  being  tested  both  experimentally  and  by 
calculation.  Using  the  second  order  Godunov  scheme,  an  idealized  setup  involving  various 
thermal  layers  had  already  been  studied  by  Glowacki,  et  al.  [32]  These  studies  clearly  in¬ 
dicated  the  possibility  of  thermally  induced  precursor  effects;  additionally,  the  calculations 
demonstrated  the  high  resolution  and  accurate  computation  of  multiple  vortex  rollups  and 
interactions  in  the  very  subsonic  region  well  behind  the  incident  shock.  This  work  was  to 
inspire  a  substantial  effort  in  the  study  of  shear  layers  during  the  contract  period  covered 
by  this  report  and  will  be  discussed  below.  In  any  event,  the  code  was  used  to  simulate 
nuclear  blast  events  via  very  long  CRAY  runs.  Some  of  this  work  is  mentioned  in  [25]  and 
references  therein  and  it  formed  the  basis  for  a  wide  variety  of  computational  simulations 
which  will  be  one  of  the  subjects  of  Section  5  below. 

Alongside  all  of  these  large-scale  computations,  a  revolution  was  occurring  in  the  CFD 
community  as  well  as  the  more  theoretically  oriented  group  involved  in  studying  schemes 
and  algorithms.  This  was  largely  driven  by  the  exponentially  increasing  power  of  the  su¬ 
percomputers,  but  the  influence  of  the  mathematical  theory  of  conservation  laws  developed 
by  Lax,  Glimm,  Liu,  DiPerna,  etc.  was  also  extremely  important.  Some  typical  examples 
of  applying  these  ideas  to  develop  schemes  may  be  found  in  the  references  [21],  [35],  [54], 
[57],  [65];  a  particularly  useful  reference  is  the  review  of  Yee  [70]  which  contains  about 
200  references  therein.  Of  course,  the  high-order,  Godunov  schemes  are  the  prototypical 
example  of  borrowing  from  the  theory  of  hyperbolic  conservation  laws.  In  any  event,  it 
is  now  the  case  that  1-D  gas  dynamics  problems  can  be  considered  solved  by  at  least  20 
different  variants  of  the  idea  of  second-  or  higher-order  upwinding. 

Additionally,  the  second-order  Godunov  technology  has  been  extended  in  many  useful 
ways  since  1985.  Examples  are  the  mesh  refinement  strategy  of  Berger  and  Colella  [4], 
the  unsplit  algorithm  of  Colella  [10]  for  multidimensional  problems  and  the  front  tracking 
method  of  Chern  and  Colella  [6],  Other  examples  involving  myself  are  discussed  below. 
At  the  present  time,  these  efforts  are  concentiated  at  LLNL  and  are  being  organized  by 
P.  Colella  and  J.  Bell. 

After  a  review  and  update  on  self-similar  oblique  shock  wave  reflection  in  Section  1, 
the  next  three  sections  discuss  progress  in  various  code  development  projects  that  I  am 
involved  with  under  this  contract.:  Section  5  lists  the  many  applications  that  I  have  been 
working  on  with  varying  degrees  of  effort  during  the  last  three  years;  since  all  of  this  work 
has  been  published  either  in  the  open  literature  or  as  technical  reports,  it  would  serve  no 
purpose  to  discuss  the  results  and  so  this  section  is  brief.  Finally,  along  with  J.  P.  Collins, 
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J.  Krispin  and  others,  I  have  involved  myself  with  various  aspects  of  isentropic  gas  dynamics 
and  testing  different  algorithms  in  2-D  settings  for  these  equations;  this  is  the  subject  of 
Section  6  which  will  also  be  brief  since  the  results  will  soon  be  available  in  preprint  form. 
Related  to  the  work  of  Section  6,  Krispin  and  myself  have  undertaken  an  analogous  study 
of  various  versions  of  the  second  order  Godunov  scheme  as  applied  to  3-D  supersonic  steady 
gas  dynamics  [40];  3-D  calculations  were  achieved  by  operator  splitting  using  the  algorithm 
of  Glaz  and  Wardlaw  [31].  Although  this  work  was  not  funded  by  DNA,  it  may  still  be  of 
interest  to  some  readers. 


I. 1  SELF-SIMILAR  SHOCK  WAVE  REFLECTION. 

Self-similar  oblique  shock  wave  reflection  for  the  equations  of  unsteady  inviscid  gas 
dynamics  is  among  the  outstanding  unsolved  problems  in  nonlinear  partial  differential 
equations  (PPE).  This  is  so  due  not  only  to  its  standing  as  the  simplest  possible  nontrivial 
shock  wave  interaction  problem  for  these  equations  and  the  importance  of  the  problem 
in  engineering  applications,  but  to  the  tantalizingly  large  amount  of  experimental  and 
computational  data  which  is  available  after  half  a  century  of  intensive  research  since  WW 

II, 

The  equations  of  unsteady,  inviscid  gasdynamics  in  Cartesian  coordinates  arc 


(la)  Pt+(pu)z+(pv)y-0 

(lb)  (pu)t  +  (pu2  +  p)x  +  (puv)y  =  0 

(lc)  (pv),  +  (puv)x  +  (pv2  +  p)y  =  0 

(ld)  ( pE)t  +  (puE  +  up)z  +  (pvE  +  vp)y  =  0 


where  p  is  the  density,  u  =  (u,v)  is  the  velocity  field,  E  =  |(u2  +  v2)  +  e  is  the  total  specific 
energy,  e  is  the  specific  internal  energy,  and  p  is  the  pressure.  The  system  is  closed  by 
specifying  an  equation-of-state  (EOS),  p  =  p (p,e);  an  important  example  is  the  polytropic 
EOS,  p  =  (7  >  l)pe,  where  7  >  1  is  constant.  Equations  (1)  are  a  system  of  conservation 
laws 

U,  +  F(U)Z  +  G(U)y  =  0 

where  U  =  (p,  pu ,  pv ,  pE)f,  etc.  The  initial  conditions  for  oblique  shock  wave  reflection  are 
set  up  as  follows:  an  incident  shock  wave,  I,  of  shock  wave  Mach  number  Ms  traverses  an 
ambient  (;.e.,  Uo  =  (0,0))  medium  and  approaches  a  wedge  surface  inclined  at  0W  degrees; 
since  the  upstream  flow  is  ambient,  M,  =  ct/'cq  where  <7  =  shock  speed  and  eg  is  the 
upstream  sound  speed.  The  coordinate  system  is  set  by  taking  the  corner  at  (x,y)  =  (0,0) 
and  the  time  t  —  0  at  the  instant  the  shock  reaches  the  corner.  There  are  no  length  scales 
in  the  initial  data,  the  equations  (1),  or  the  EOS,  so  it  is  reasonable  to  look  for  self-similar 
or  pseudosteady  solutions  for  t  >  0  in  the  similarity  variables 

(2)  (£,»j)  -  (x/t,y/t).. 

In  conservation  form,  the  transformed  equations  are 


(3a)  (pu)i  +  (pv)„  =  -2  p 

(3b)  (pu2  +  p){  +  ( puv)v  =  -3  pu 

(3c)  (puv)t  +  (pv2  -f  p)v  =  -3 pv 

(3d)  (puH)^  +  (pvH)v  =  —p(u2  -f  v2)  —  2pH 
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where  Ci  =  (u,  v)  is  the  self-similar  velocity  field,  u  =  u  —  £,  v  =  v  —  r?,  H  —  u 2  +  v2)  +  h 

is  the  self-similar  total  enthalpy,  and  h  =  e  +  p/p  is  the  specific  enthalpy.  Additionally, 

(4)  M2  —  (u2  +  v2)/c2 

defines  the  self-similar  Mach  number  M,  where  c  =  sound  speed.  In  vector  form,  the 
equations  (3)  are  written 

+  Gjf  =  § 

where  F  —  ( pu ,  pu 2  +p,  put),  puH )*,  etc.  One  notes  that  the  first-order  part  of  this  system, 
4  (5,  =  0,  arc  exactly  the  equations  of  steady  gas  dynamics  in  the  (£,»?)  coordinates. 

The  wall  boundary  condition  for  equations  (1)  is  perfect  reflection,  i.e.,  flow  is  parallel 
to  the  wall.  It  is  easily  seen  that  this  condition  transforms  analogously  for  equations 
(3).;  The  boundary  condition  at  infinity  for  equations  (1)  is  the  time-dependent  Dirichlet 
condition  provided  by  the  Rankine-Hugoniot  conditions  connecting  states  U0  and  U\.  Since 
no  disturbance  can  reach  infinity  in  finite  time,  this  boundary  condition  can  also  be  easily 
transformed.  Further,  the  boundary  may  be  moved  to  a  finite  distance  from  the  origin 
so  long  as  it  is  outside  the  region  of  disturbed  flow;  the  intersection  of  the  incident  shock 
with  such  a  boundary  is  easily  calculated  in  the  (£,  r;) -frame  for  t  >  0  by  using  the  known 
shock  speed,  a. 

Unlike  equations  (1),  which  are  always  hyperbolic  under  very  general  conditions  on  the 
EOS,  the  equations  (3)  of  steady  gas  dynamics  with  source  terms  are  hyperbolic  only  when 
the  flow  is  supersonic  in  the  (£,  T])- frame.  For  subsonic  flow,  the  streamline  characteristic 
remains  real  but  the  sound  wave  characteristics  become  complex.  Thus,  equations  (3)  are 
of  mixed  type.  Note,  however,  that  for  the  boundary  conditions  under  consideration  here, 
the  flowfield  becomes  supersonic  as  (£,  rj)  —*  oo  for  any  initial  data;  in  particular,  any  finite 
farfield  boundary  chosen  as  above  should  satisfy  the  condition  that  the  flow  is  supersonic 
nearby. 

Summarizing,  we  are  interested  in  solving  the  boundary  value  problem  given  by  equa¬ 
tions  (3)  subject  to  the  given  boundary  conditions.  Evidently,  the  solution  depends  para¬ 
metrically  on  Ma,0w,  the  EOS  and  the  ambient  data  (po,eo).  For  the  important  case 
of  a  polytropic  gas,  the  parameters  are  MB,8W  and  7;  the  ambient  state  is  no  longer  an 
independent  parameter.  To  the  best  knowledge  of  the  author,  there  do  not  exist  any  math¬ 
ematical  results  pertaining  to  the  question  of  globed  solutions  for  this  problem;  the  issues 
of  existence,  uniqueness,  and  continuous  dependence  on  the  parameters  remain  open  for 
all  parameter  l  unges. 

Any  conjectures  that  we  have  concerning  these  questions  must  be  inferred  from  shock 
tube  experiments  and  CFD  simulations.  Additionally,  mathematical  analysis  is  often  ap¬ 
plied  to  infinitesimal  or  local  flowfield  regions,  especially  shock  wave  interactions,  once 
such  regions  are  apparent  from  experiments  or  calculations.  Of  course,  such  experiments 
and  calculations  are  subject  to  experimental  error  (  e.g.,  the  incident  shock  has  a  finite 
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thickness  and  may  not  move  at  precisely  constant  speed,  the  wall  will  not  be  perfectly 
smooth,  data  analysis  errors,  etc.)  and  numerical  error.  Abo,  any  experimental  result 
will  represent  a  solution,  not  of  equations  (1),  but  of  an  augmented  version  of  equations 
(1)  which  takes  into  account  real  gas  effects  such  as  vibrational  relaxation,  dissociation- 
recombination  chemistry,  and  viscosity.  In  some  examples  (especially  at  high  Ma),  these 
effects  are  significant;  in  all  cases,  it  is  necessary  to  carefully  analyze  them.  Numerical 
simulations  may  not  be  fully  resolved,  may  introduce  ‘artificial’  viscosity,  and  may  contain 
other  systematic  errors;  furthermore,  it  is  not  possible  to  prove  convergence  or  obtain  use¬ 
ful  error  estimates  at  the  present  time.  Despite  all  of  these  caveats,  most  workers  in  the 
field  would  at  least  maintain  that  the  data  is  very  useful  in  deducing  properties  of  solutions 
of  equations  (3).  Still,  it  must  be  emphasized  that  the  answers  to  these  questions  are  not 
at  all  obvious. 

For  further  analytical/mathematical  analysis  of  self-similar  oblique  shock  wave  reflec¬ 
tion,  the  reader  may  find  [26]  useful.  Recently,  several  review  articles  on  this  subject  have 
appeared;  we  mention  those  by  G.  Ben-Dor  [3],  1. 1.  Glass  [23]  and  H.  Hornung  [39].,  These 
references  also  lead  to  many  hundreds  of  papers  in  the  field.  We  single  out  the  very  im¬ 
portant  work  of  L.  F.  Henderson,  et  al.  [36],  [37],  [38]  on  the  subject  of  transition  between 
regular  and  Mach  reflection  which  is  the  outstanding  open  question. 

As  noted  in  the  Introduction,  this  problem  has  been  the  subject  of  many  numerical 
studies.;  This  is  due  not  only  to  it’s  inherent  interest,  but  also  to  it’s  suitability  as  a 
difficult  test  problem  for  numerical  methods.  The  earliest  work  of  this  type  dates  back  to 
the  mid-1970’s  and  was  done  in  the  aerospace  CFD  community;  see  [51],  [58]  and  [59]. 

An  extension  of  the  basic  second  order  Godunov  technology  was  developed  for  the  high 
Ms  cases  where  real  gas  effects  in  low  density  air  lead  one  to  suspect  that  the  flowfield  is  not 
in  equilibrium  (and,  therefore,  is  not  self-similar  since  the  relaxation  times  introduce  new 
scales  into  the  equations).  The  results  were  finalized  and  appeared  in  the  open  literature 
during  the  contract  period,  see  Glaz,  et  al.  [27].  These  new  calculations  were  all  that 
we  could  have  hoped  for  since  they  explained  away  all  of  the  quantitative  disagreements 
discussed  in  the  precedent  work  [28],  [29]  where  nonequilibrium  effects  were  conjectured 
to  be  decisive.  It  may  also  be  mentioned  that  a  code  of  this  type,  adapted  for  general 
geometries,  would  be  quite  useful  in  high  altitude  aerospace  studies. 

Finally,  a  very  nice  numerical  study  of  the  weak  limit  case  has  been  undertaken  by 
Colella  and  Henderson  [14].  The  numerical  results  therein  seem  to  have  resolved  some  of 
the  open  questions  among  the  large  subcommunity  interested  in  this  problem. 
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1.2  A  NEW  PROJECTION  METHOD  FOR  INCOMPRESSIBLE  FLOW. 

Around  1968,  A.  J.  Chorin  and  R.  Temam  developed  projection  methods  for  the  in¬ 
compressible  Euler  and  Navier-Stokes  equations.  This  is  an  inherently  split  method  in 
which  the  advection  or  advection-diffusion  step  takes  the  velocity  field  out  of  the  space 
of  divergence-free  vector  fields  (this  may  be  thought  of  either  at  the  PDE  level  or  at  the 
discrete  level  with  suitable  discrete  divergence  and  gradient  operators)  and  the  projection 
step  calculates  the  divergence-free  part  of  this  field.  Of  course,  this  is  a  linear  algebra 
problem.  Bell,  et  al.  (1]  have  updated  this  algorithm  in  three  ways:  (1)  the  linear  algebra 
has  been  significantly  modernized,  (2)  the  coupling  between  the  two  split  steps  is  stronger 
and  (3)  the  advection  algorithm  has  been  radically  altered  from  a  standard  finite  differ¬ 
ence  to  a  version  of  the  unsplit  second  order  Godunov  scheme.  The  latter  step  has  led  to 
spectacular  improvements  in  the  resolution  capabilities  of  the  method. 

Further  work  in  extending  the  algorithm  for  the  Boussinesq  approximation,  fully  vari¬ 
able  density  flow  and  3-D  versions  is  continuing  at  LLNL  by  J.  Bell  and  colleagues.  Ulti¬ 
mately,  the  method  may  be  able  to  handle  the  low-speed  part  of  fireball  calculations  and 
couple  with  a  compressible  algorithm  at  some  point  for  the  late  stages.  In  this  sense,  it 
is  a  competitor  for  the  implicit-explicit  technology  discussed  in  Section  3.  Also,  it  has 
occurred  to  many  of  us  that  the  zero  Mach  number  limit  of  the  latter  algorithm  may  well 
be  a  version  of  the  new  projection  algorithm  and  that  understanding  this  limit  may  be  the 
key  to  obtaining  sharp  estimates  on  the  performance  of  both  schenv  classes. 
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1.3  IMPLICIT-EXPLICIT  SCHEMES. 

A  hybrid  implicit-explicit  scheme  of  a  type  first  discussed  in  [22]  for  the  case  of  La- 
grangian  hydrodynamics  has  been  developed  for  the  equations  of  Eulerian  hydrodynamics 
by  Collins,  et  al.  [18].  For  these  schemes,  the  difference  approximation  in  time  is  either 
implicit  or  explicit,  separately  for  each  family  of  characteristics  and  for  each  cell  in  the  fi¬ 
nite  difference  grid,  depending  on  whether  the  local  CFL  number  for  that  family  is  greater 
than  or  less  than  one..  In  both  cases  the  hybridization  is  continuous  at  CFL  number  equal 
to  one,  and  the  scheme  for  the  explicit  modes  is  a  second-order  Godunov  method  of  a  type 
discussed  in  [16],  [13]. 

This  implicit-explicit  strategy  is  intended  for  problems  with  spatially  and/or  tempo¬ 
rally  localized  stiffness  in  wave  speeds.  By  stiffness,  we  mean  that  the  high  speed  modes 
contain  very  little  energy,  yet  they  determine  the  explicit  time  step  through  the  CFL  con¬ 
dition,  For  hydrodynamics,  the  main  example  is  nearly  incompressible  flow.  Here,  the 
sonic  waves  will  be  nearly  acoustic  and  largely  decoupled  from  the  particle  modes.  Tra¬ 
ditionally,  such  problems  are  handled  in  one  of  two  ways:,  at  the  PDE  level  by  solving 
instead  the  incompressible  limit  of  the  set  of  governing  equations,  or  numerically  by  artifi¬ 
cially  increasing  the  temperature  of  the  gas  in  such  a  way  that  the  Mach  number  becomes 
significant  but  still  low  enough  that  compressibility  effects  are  negligible. 

Another  example  is  provided  by  magnetohydrodynamic  calculations  in  certain  config¬ 
urations  [5],  [56].  Here,  if  the  problem  leads  to  spatially  localized  regions  where  a  low 
density  plasma  is  immersed  in  a  strong  magnetic  field  then  the  Alfven  wave  speed  can  be¬ 
come  very  large.  Our  approach  can  potentially  deal  with  this  problem  in  a  straightforward 
manner  and  avoid  ad  hoc  solutions,  e.g.,  modifying  the  calculation  of  the  electric  field  at 
low  densities. 

A  related  example  is  the  problem  of  shock  wave  -  boundary  layer  interaction.  It  is 
often  important  to  simultaneously  perform  a  high  resolution  calculation  of  the  strong  wave 
interactions  in  the  free  stream  while  also  resolving  enough  of  the  detail  in  the  boundary 
layer  to  accurately  model  the  effect  of  the  no-slip  boundary  condition  on  the  free  stream 
as  well  as  handling  any  transport  effects  between  the  boundary  and  the  inviscid  flow  region 
(which  might  include  particle  and  temperature  transport  from  the  boundary  condition). 
The  boundary  layer  length  scale  is  usually  several  orders  of  magnitude  smaller  than  even 
a  single  computational  zone  width,  while  the  particle  speed  is  approaching  zero  along  with 
the  strengths  of  the  sonic  waves  as  the  boundary  is  approached.  These  difficulties  can  be 
partially  ameliorated  by  constructing  appropriate  adaptive  meshes,  but  this  is  not  usually 
sufficient.;  A  common  approach,  especially  in  engineering  calculations  where  efficiency  is 
important,  is  not  to  use  boundary  layer  zoning  at  all  but  to  explicitly  difference  some 
approximation  to  the  Navier-Stokes  equations  in  the  first  several  zones  near  the  boundary, 
thereby  satisfying  the  correct  boundary  condition.  Of  course,  accuracy  can  be  a  problem 
in  such  calculations.  Alternatively,  one  can  use  boundary  layer  zoning  at  the  expense  of 
also  differencing  the  entire  set  of  equations  implicitly,  at  least  in  a  locally  refined  region, 
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in  order  to  avoid  stability  problems;  a  discussion  of  such  methods  for  upwind  schemes  may 
be  found  in  [70].  Potentially,  the  latter  approach  is  viable,  but  the  interface  between  the 
boundary  layer  and  inviscid  regions  can  be  a  difficult  numerical  problem.  The  extension  of 
our  approach  to  the  Navier-Stokes  equations,  also  developed  in  [18],  allows  for  the  smooth 
transition  between  a  high  resolution  explicit  scheme  in  the  exterior  flowfield  and  a  stable 
method  in  the  boundary  layer. 

The  development  of  our  implicit-explicit  scheme  has  followed  the  following  design 
principles:; 

(1)  The  scheme  is  in  conservation  form,  both  for  conservation  laws  and  for  their  viscous 
extensions. 

(2)  The  implicit-explicit  hybridisation  is  a  continuous  switch  and  operates  on  each 
characteristic  field  separately;  the  method  is  entirely  local  in  both  time  and  space, 
depending  only  on  the  data  in  each  computational  zone. 

(3)  In  the  event  that  all  characteristic  modes  are  explicit,  the  scheme  reduces  to  a 
version  of  the  second  order  Godunov  scheme  [16],  [13]. 

(4)  The  implicit  advection  scheme,  which  is  only  first  order  accurate  in  time,  satisfies 
a  maximum  principle  and  is  unconditionally  stable. 

(5)  In  the  limit  of  steady  state,  the  scheme  is  second  order  accurate  in  space. 

The  first  three  of  these  properties  are  satisfied  by  the  scheme  of  Fryxell,  et  al.  [22].  Their 
implicit  scheme  is  also  second  order  accurate  in  time,  but  does  not  satisfy  a  maximum 
principle  (and,  therefore,  is  not  monotone).;  Additionally,  for  a  system  of  n  equations,  their 
method  requires  the  inversion  of  a  block  2n  x  2n  tridiagonal  system  whereas  the  current 
method  requires  only  a  block  n  x  n  tridiagonal  system.  For  our  intended  applications,  we 
do  not  regard  the  lack  of  second  order  temporal  accuracy  to  be  a  serious  problem.  Waves 
that  we  wish  to  resolve  in  a  flowfield  will  be  treated  using  the  second  order  accurate  explicit 
scheme;  the  others  will  be  rapidly  relaxed  to  equilibrium. 

A  key  component  of  the  work  described  in  [18]  is  the  introduction  of  an  appropriately 
smooth  numerical  flux  function  for  the  hyperbolic  equations.  We  use  a  suitable  version 
of  the  Engquist-Osher  flux  [21]  as  modified  for  systems  by  Bell,  et  al.  [2].  The  version 
used  is  sufficiently  smooth  so  that  the  Newton’s  method  linearization  is  well-behaved;  in 
particular,  it  converges  to  steady  states  even  in  the  presence  of  strong  shocks. 

The  special  case  of  one-dimensional  inviscid  and  viscous  compressible  flow  is  treated 
in  [18]  and  a  polytropic  equation-of-state  is  assumed  for  convenience.  However,  the  ideas 
and  implementation  methodology  are  easily  extended  to  more  general  systems  in  one  space 
dimension.  In  particular,  it  is  expected  that  our  results  can  be  extended  to  gas  dynamics 
with  a  general  equation-of-state  using  the  ideas  in  [13].  A  more  difficult  question  is  the 
extension  of  these  ideas  to  two  or  three  space  dimensions.  A  few  operator  split  2D  boundary 
layer  calculations  have  already  been  successfully  completed  by  Collins  [17]. 

At  the  present  time,  Collins  is  attempting  to  develop  an  unsplit  version  of  the  method- 
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ology  using  the  explicit  unsplit,  second-order  Godunov  scheme  introduced  by  Colella  [10] 
as  a  starting  point;  this  work,  if  successful,  will  lead  to  his  Ph.  D.  degree  in  Applied 
Mathematics  at  the  University  of  Maryland,  College  Park  (I  am  his  thesis  advisor).  Our 
strategy  is  to  concentrate  first  on  the  zero  Mach  number  limit  for  gas  dynamics  with  peri¬ 
odic  boundary  conditions;  as  noted  in  Section  2,  this  approach  may  well  lead  to  significant 
new  ideas  in  several  directions.  In  addition,  we  are  experimenting  with  alternatives  to 
Newton  linearization  such  as  multigrid  acceleration  and  preconditioned  conjugate  gradi¬ 
ents  applied  to  the  discrete  nonlinear  problem  for  the  implicit  modes.  Later,  nontrivial 
boundary  conditions  will  be  introduced  along  with  various  approximations  of  the  viscous 
terms  in  the  Navier-Stokes  equations  so  that  boundary  layer  transport  can  be  calculated. 
Ultimately,  some  version  of  mesh  refinement  will  be  required  so  that  the  full  Navier-Stokes 
equations  can  be  solved  in  critical  regions,  e.g.,  the  reflection  point  in  studies  of  transition 
between  regular  and  Mach  reflection.  When  this  point  is  reached,  it  is  fortunate  that 
experimental  work  from  UTIAS  will  be  available  which  may  well  be  useful  in  validating 
and/or  improving  the  algorithm;  we  refer  to  the  papers  [61],  [62],  [66]  in  which  transition 
was  studied  as  a  function  of  Reynolds  number. 
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1.4  MULTIMATERIAL  ALGORITHMS. 


Colella,  et  al.  [11],  have  developed  a  volume-of-fluid  type  method  for  the  numerical 
calculation  of  compressible  flow  problems  in  which  the  fluid  is  made  up  of  a  number  of 
thermodynamically  distinct  materials  separated  by  sharp  interfaces.  In  this  approach,  a 
standard  finite  difference  representation  of  the  solution  is  augmented  by  cell-centered  val¬ 
ues  for  the  thermodynamic  quantities:  pa ,  e“,  Va,  are  the  density,  internal  energy,  and 
fractional  volume  occupied  by  the  atk  fluid  in  each  zone,  a  =  1 ,  ,  n/.  The  evolution 
of  this  representation  can  be  thought  of  as  consisting  of  two  parts.  One  is  the  effective 
Lagrangian  dynamics — the  accelerations,  compressions,  and  work  done  on  the  multifluid 
representation — which  is  computed  under  the  assumption  that  the  various  fluid  compo¬ 
nents  are  in  pressure  equilibrium  with  one  another  in  each  cell,  and  that  each  cell  has 
a  single  velocity.  From  a  physical  point  of  view,  the  assumption  of  pressure  equilibrium 
is  not  unreasonable,  since  the  pressure  is  continuous  across  a  contact  discontinuity.  The 
requirement  that  the  cell  has  a  single  velocity  is  not  an  appropriate  one  in  more  than  one 
dimension,  since  slip  can  be  generated  at  a  fluid  interface.  Thus,  the  jump  in  the  ther¬ 
modynamic  variables  across  the  interface  is  tracked,  while  the  jump  in  tangential  velocity 
is  captured  using  the  underlying  conservative  finite  difference  method.  The  other  part  of 
the  evolution  is  the  motion  of  the  fluid-interface  through  the  finite  difference  grid.  This 
is  done  by  reconstructing  locally  the  interface  geometry  l'rom  the  fractional  volumes,  and 
transplanting  material  along  streamlines  defined  by  the  single  fluid  velocities..  Multifluid 
methods  have  been  in  use  for  some  time  as  noted  in  the  Introduction  and  have  been  quite 
effective  in  representing  complicated  multifluid  configurations  undergoing  large  distortions. 

Two  innovations  into  this  class  of  algorithms  are  introduced  in  [11]..  The  first  one 
concerns  the  effective  Lagrangian  dynamics  of  the  multifluid  cells.  A  formulation  of  this 
dynamics  which  is  thermodynamically  consistent  is  derived  there.  If  the  various  fluid  com¬ 
ponents  in  a  cell  are  in  pressure  equilibrium,  then  they  remain  so  to  leading  order  in  the 
truncation  error,  assuming  the  the  pressure  gradients  and  compressions  art  finite.  Since 
that  assumption  will  occasionally  be  violated-for  example,  when  a  shock  over-takes  a  ma¬ 
terial  interface-a  relaxation  scheme  to  restore  pressure  equilibrium  in  multifluid  cells  is 
introduced.  The  second  innovation  is  the  coupling  of  this  method  to  an  operator-split, 
second-order,  Eulerian,  Godunov  method.  In  previous  multifluid  algorithms,  the  concep¬ 
tual  division  into  two  parts  was  used  literally  as  the  basis  for  the  design  of  the  algorithm, 
with  the  underlying  single-fluid  algorithm  having  to  be  formulated  as  a  Lagrangian  step 
followed  by  a  remap.  In  the  approach  taken  in  [11],  the  underlying  difference  algorithm  is 
a  conservative  predictor-corrector  method  which  provides  pressures  and  velocities  at  the 
cell  edges  as  part  of  the  flux  calculation. 

Both  1-D  and  2-D  calculations  are  presented  in  [11].  In  fact,  2-D  supersonic  jets  have 
been  very  successfully  calculated  with  previous  versions  of  the  algorithm  (in  a  few  cases, 
the  results  indicated  new  physics  not  previously  calculated  by  others  working  the  same 
problem  -  a  set  of  models  based  on  astrophysical  jets).  A  key  change  in  the  new  version  is 
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a  slight  modification  in  the  technique  used  by  the  underlying  single  fluid  algorithm  in  it’s 
handling  of  general  EOS.  This  change  has  enabled  us  to  perform  HE  product  gas/water 
calculations  using  the  Gittings  EOS  for  water.  It  is  hoped  that  the  new  scheme  will  find 
applications  in  this  environment  for  the  DNA  community. 

Finally,  a  very  successful  study  has  been  undertaken  by  Colella  et  al  [15]  in  simulating 
several  cases  of  shock  wave  refraction  and,  in  some  runs,  comparing  with  experimental 
work  of  L.  F.  Henderson.  These  calculations  were  based  on  a  version  of  the  algorithm 
presented  in  [11]. : 
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1.5  APPLICATIONS. 

A  major  effort  has  been  undertaken  over  the  past  hve  years  in  studying  free  shear  layers 
in  an  attempt  to  prove  that  instability  and  rollup  is  a  purely  inviscid  dynamic  phenomenon 
and  is  independent  of  the  Reynolds  number,  see  [7-9]  and,  especially,  [43].  I  believe  that 
our  simulation  of  the  famous  Brown-Roshko  experiment  of  1974  which  is  reported  in  great 
detail  in  [43]  is  the  best  calculation  ever  made  using  the  old  split  version  of  the  second 
order  Godunov  scheme.  It  is  clear  that  major  contributions  have  been  made  to  this  area 
of  fluid  mechanics. 

The  reports  [42],  [44-50]  detail  the  many  calculations  performed  from  NSWC  over  the 
past  few  years  in  studying  wall  oounded  shear  layers  and  jets,  precursors  and  decursors  for 
both  simulations  and  airblast  modeling  and  gas-particle  2-phase  boundary  layer  modeling 
under  the  restraints  involved  in  using  an  inviscid  code.  Dense  gas  models  were  used  as  well 
as  the  dusty  gas  limit  model,  see  [53].  Many  of  these  calculations  were  very  large-scale, 
involving  massive  data  sets  and  complex  postprocessing  and  analysis. 

Finally,  I  was  involved  in  a  demonstration  study  using  the  new  quadrilateral  unsplit 
code  developed  at  LLNL;  we  set  up  various  loading  configurations  already  studied  ex¬ 
perimentally  at  UTIAS  and  male  several  very  successful  comparison  calculations  [24]. 
Currently,  we  are  now  redoing  these  runs  using  the  even  newer  quadrilateral  AMR  code; 
the  preliminary  results  are  quite  spectacular  and  we  would  like  to  finish  this  study  soon. 
The  reader  may  be  interested  in  earlier  loading  calculations  using  the  operator  split  version 
of  the  code;  see  [25]  and  references  therein. 
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1.6  ISENTROPIC  GAS  DYNAMICS. 

The  equations  of  2D  unsteady  isentropic  gas  dynamics  are  given  by  equations  (la), 
(lb),  and  (1c)  of  Section  1  with  the  EOS  given  by  the  relation  p  =  p(p).  Just  as  for  gas 
dynamics,  this  set  of  equations  forms  a  system  of  hyperbolic  conservation  laws 

Ut  +  F(U)xj-G(U)y  =  0 

where  U  =  (p,  pu,  pv)1,  etc.  This  effectively  replaces  the  energy  equation  by  the  constraint 
that  the  entropy  be  conserved  along  streamlines;  in  practice,  we  have  assumed  that  the 
entropy  is  everywhere  a  constant.  These  equations  are  valid  in  two  limiting  situations. 
First,  since  the  jump  in  entropy  across  a  discontinuity  is  of  third  order  in  the  shock  strength 
and  entropy  is  conserved  in  smooth  flow,  this  set  is  an  excellent  approximation  if  all  shocks 
may  be  assumed  to  be  weak  a  priori.  Also,  the  system  is  equivalent  to  a  version  of  nonlinear 
acoustics  if  the  variables  are  reinterpreted  appropriately.  The  same  is  true  for  the  shallow 
water  equations.  Second,  the  entropy  jump  is  small,  even  for  very  strong  waves,  if  the 
material  density  is  large  and  the  EOS  is  stiff;  the  example  of  inijrest  here  is  shock  loading 
of  metals  although  some  of  the  ideas  may  also  be  useful  for  water. 

In  one  space  dimension  the  fluid  dynamic  equations  are  equivalent  to  the  equations 
of  solid  mechanics  in  their  simplest  formulation  with  the  exception  that  the  pressure  is 
viewed  as  the  stress  but  of  opposite  sign;  of  course,  the  constitutive  modeling  (i.e.,  the 
EOS  here)  may  be  more  difficult.  Some  materials,  such  as  iron  and  Plexiglas,  undergo  a 
phase  change  when  shocked  to  a  sufficiently  high  pressure.  To  model  this  phenomenon, 
it  is  necessary  to  allow  the  EOS  to  be  non-convex  in  the  appropriate  pressure  ranges;  the 
theory  of  hyperbolic  conservation  laws  predicts  more  complicated  wave  structures  in  this 
case.  This  theory,  in  particular  as  it  applies  to  the  Riemann  problem,  is  well  understood 
and  agrees  with  results  from  experiment  for  those  few  situations  where  such  experiments 
are  possible.  A  major  objective  of  our  work  was  to  simulate  the  a  —  e  phase  transition  in 
Iron  which  occurs  at  about  137  kilobars  and  has  been  extensively  studied.  Our  model  is  too 
simple  to  allow  for  the  weak  elastic  precursor  wave  which  is  observed  in  some  experiments 
but  this  has  not  been  a  major  factor;  an  interesting  point  is  that  this  precursor  wave  is 
much  faster  than  the  following  main  wave  pattern  and,  if  we  had  included  it  in  our  model, 
the  implicit-explicit  ideas  of  Section  3  would  be  useful. 

Another  objective  of  this  work  was  to  test  several  variants  of  the  second-order  Go¬ 
dunov  scheme  in  the  relatively  simple  setting  of  isentropic  gas  dynamics,  but  for  materials 
where  the  solution  to  Riemann  problems  would  be  complex  and  very  expensive  to  compute 
exactly..  This  comparison  study  was  based  on  three  approaches:  the  original  algorithm  of 
Colella  and  Glaz  [13]  suitably  adapted  for  the  EOS,  the  very  general  methodology  of  Bell, 
et  al.  [2]  (which  is  referred  to  below  as  the  BCT  method)  and  an  extremely  simple  ver¬ 
sion  due  to  Davis  [20].  We  worked  with  a  model  EOS  as  well  as  the  iron  EOS  from  the 
literature  and  ID  and  2D  problems  were  studied.  The  2D  calculations  were  the  analogue 
of  the  oblique  self-similar  shock  wave  reflections  discussed  in  Section  1  for  the  equations 
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of  isentropic  gas  dynamics,  although  here  we  have  no  experimental  data  against  which  to 
compare. 

The  results  of  this  work  are  detailed  in  the  paper  of  Krispin,  et  al.  [41].  A  very 
brief  summary  follows.  First,  the  new  phenomenology  resulting  from  the  nonconvexity  is 
quite  exciting  and  it  would  be  of  interest  to  see  if  similar  results  are  in  the  experimental 
literature.  FVom  the  numerical  point  of  view,  the  original  algorithm  worked  quite  well  and 
was  very  expensive,  as  expected.  The  Davis  approach  is  reasonable  for  simple  solutions 
such  as  regular  reflection  but  appears  to  break  down  when  there  are  strong  2-D  wave 
interactions,  e.g.,  double  Mach  reflection;  on  the  other  hand,  it  is  a  very  efficient  method. 
The  main  conclusion  was  that  several  variants  of  the  BCT  method  yielded  excellent  quality 
results  at  an  intermediate,  but  acceptable,  expense.  We  believe  that  this  type  of  study  will 
become  more  common  in  the  field  as  shock  waves  in  exotic  materials  and  more  complicated 
systems  of  conservation  laws  are  modeled  numerically. 

Another  study  has  also  been  carried  out  with  the  same  code  for  the  simpler  case  of  a 
convex  EOS.  Here,  we  were  working  from  a  set  of  conjectures  concerning  possible  solutions 
of  the  two-dimensional  Riemann  problem  for  isentropic  gas  dynamics.  The  initial  data  was 
chosen  from  the  special  cases  in  which  each  of  the  four  one-dimensional  Riemann  problems 
has  a  solution  consisting  of  one  wave  only;  still,  the  analytic  problem  is  quite  formidable. 
The  general  case  is  surely  impossible  to  solve  other  than  by  numerical  means  so  our  efforts 
here  are  partially  a  validation  study  for  a  future  attack  on  the  general  case.  The  results 
have  been  collated  and  analyzed  with  respect  to  the  a  priori  conjectures,  see  Collins,  et  al, 
[19],  They  may  be  stated  quite  simply:  (1)  some  of  the  conjectures  are  true  and  others 
are  false  and  (2)  many  of  the  wave  interaction  patterns  are  quite  striking.  Indeed,  it  is 
possible  to  recover  all  of  the  self-similar  wave  patterns  familiar  from  the  studies  mentioned 
in  Section  1  through  judicious  choice  of  initial  data. 
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SECTION  2 
CONCLUSIONS 

I  believe  that  the  following  conclusions  are  quite  clear: 

(1)  The  older  operator  split  code  for  gas  dynamics  with  a  general  convex  equation 
of  state  along  with  the  many  driver  programs  written  to  use  it  (e.g.,  blast  wave, 
wedge,  simulator,  shear  layer,  etc.  drivers)  is  extremely  robust  and,  if  not  pushed 
too  far,  reliable.  The  results  obtained,  both  for  theoretical  physics  studies  as  well 
as  for  engineering  applications,  are  clearly  the  equal  or  better  than  those  obtained 
with  competing  CFD  technologies. 

(2)  The  more  recently  developed  versions  which  include  mesh  refinement  capability, 
general  quadrilateral  meshes,  etc.  have  extended  the  range  of  problems  which  we 
can  study  as  well  as  improve  our  studies  for  problems  for  which  both  the  older  and 
newer  methods  are  applicable  in  the  sense  that  higher  resolution  can  be  obtained 
more  efficiently  (however,  this  is  not  always  the  case), 

(3)  The  implicit-explicit  methodology  has  passed  it’s  initial  testing  and  it  seems  likely 
that  a  general  capability  will  be  available  in  a  few  years  for  several  types  of  problems 
of  great  interest  to  the  DNA  community.  The  relationship  between  this  idea  and 
the  new  version  of  the  projection  code  for  incompressible  flow  may  prove  to  be  the 
theoretical  key  in  making  fast  progress, 

(4)  The  multimaterial  extension  of  the  second  order  Godunov  scheme  is  now  available 
for  2-D  problems,  at  least  for  gases  where  it  has  been  validated  against  experimental 
data  and,  in  a  different  set  of  computations,  has  compared  favorably  with  other 
numerical  studies.  The  more  difficult  case  of  gas-water  interfaces  is  decidedly  more 
promising  than  even  a  year  ago  and  it  appears  that  2-D  applications  will  be  feasible 
with  a  little  additional  work. 

(5)  Our  studies  on  isentropic  gas  dynamics  have  indicated  that  a  great  deal  of  funda¬ 
mental  work  remains  to  be  done  in  the  area  of  numerical  schemes  for  situations 
involving  unusual  wave  pattern  phenomenology.  These  problems  -  e.g.,  solid  me¬ 
chanics,  multiphase  flow,  MHD,  etc.  -  are  coming  under  increasing  scrutiny  in  the 
mathematical  community  and  it  seems  that  the  DNA  community  has  substantial 
applications  in  these  fields.  It  should  also  be  emphasized  that  the  phenomenology 
is  of  interest  in  its  own  right. 
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SECTION  3 

RECOMMENDATIONS 


(1)  A  greater  effort  should  be  made  to  make  proven  codes  such  as  the  original  op¬ 
erator  split,  second-order  Godunov  code  for  general  convex  EOS  available  to  the 
community  at  large.  I  do  not  recommend  making  drivers  freely  available  except  on 
a  case-by-case  basis.  The  same  goes  for  new  ard/or  experimental  codes. 

(2)  I  strongly  support  DNA’s  efforts  at  the  central  site  in  the  area  of  graphics  support. 

(3)  It  seems  that  personal  workstations  are  becoming  more  powerful  at  an  accelerating 
rate;  DNA  applications  may  soon  be  possible  on  such  equipment.  On  the  other 
hand,  it  may  become  increasingly  possible  to  perform  3-D  computations  on  mas¬ 
sively  parallel  machines  (although  I  believe  that  this  possibility  is  being  overrated 
in  the  near  term,  at  least  in  terms  of  solving  significant  problems  such  as  transition 
and  secondary  transition  to  3-D  turbulence  followed  by  highly  resolved  simulations 
of  the  resulting  fully  turbulent  flowfield).  Thus,  it  is  at  least  possible  that  CRAY 
equipment  may  be  squeezed  from  both  ends.  This  should  be  studied  in  about  two 
years  and  resolved  soon  thereafter  from  the  point  of  view  of  DN  A  interests. 

(4)  The  very  substantial  code  development  efforts,  centered  at  LLNL,  have  paid  off  in 
many  important  ways  and  should  be  continued.  I  believe  that  the  further  study  of 
the  implicit-explicit  idea  is  highly  warranted  in  view  of  the  large  payoff  if  it  works 
out. 

(5)  Somebody  should  find  out  whether  or  not  second-order  Godunov  schemes  work 
on  triangular  meshes  (by  which  I  mean  general  and,  preferably,  moving  meshes, 
and  not  simply  subdivisions  of  rectangular  meshes  for  which  special  treatments  are 
clearly  possible). 

(6)  Shock  waves  in  multiphase  flow  is  a  very  interesting  mathematical  problem  as 
evidenced  by  the  work  of  Wendroff  [66-67]  and  see  also  [55].  A  study  should  be 
made  to  determine  whether  or  not  numerical  scheme  development  is  warranted  in 
this  area  relative  to  DNA  interests.  On  the  other  hand,  it  is  obvious  that  codes 
are  required  for  multiphase  boundary  layers  and  I  very  strongly  recommend  that 
our  efforts  aimed  at  this  goal  be  continued. 
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